# -*- coding: utf-8 -*-
"""
Created on Mon Jan  7 10:51:22 2013
This module 
@author: daniel
"""

import numpy as np
from scipy.optimize import minimize
from scipy.optimize import minimize_scalar
from miscTools import frange
from FlightConditions import ISAtmosphere
from scipy import interpolate as interp
import time

class thrustAnalysis():
    """
    Provides functions for calculating total aircraft thrust and for finding
    engine power settings and fuel consumption such that a required thrust
    is matched. Wraps the thrust analysis modules so they can be accessed via
    aircraft objects. This class is available from aircraft.thrustModule
    
    Parameters
    ----------
    
    aircraft : object
        aircraft object. Class automatically initializes by propeller and 
        engine of the aircraft.
    
    
    Examples
    --------
    
    initialize thrust analysis module
    
    >>> import aircraft
    >>> ac = aircraft.load('sampleInput1')
    >>> th = thrustAnalysis(ac)
    """
    def __init__(self,aircraft):
        self.aircraft       =aircraft
        self.betaRange      =self.aircraft.propeller.betaRange
        self.propellerList  =self.aircraft.propeller.propellerList
        self.engineList     =self.aircraft.engine.engineList
        self.mode           ='max_thrust'
        self.nEng           =len(self.engineList)
        self.thrust         =0.0
        self.power          =0.0
        self.fuelFlow       =0.0
        self.RPM            =np.array([])
        self.propEfficiency =np.array([])
        self.beta           =np.array([])
        self.powerSetting   =np.array([])
    def runAnalysis(self,V,rho,powerSetting):
        """
        Determines the thrust of variable/fixed/single/multi engine configs
        
        Parameters
        ----------
        
        V : float, m/sec
            velocity
        rho : float, kg/m**3
            air density
        powerSetting : float [0;1]
            power setting
        
        Examples
        --------
        
        Single engine variable pitch
    
        >>> velocity = 60.
        >>> density = 1.2255
        >>> pSetting = 100.0
        >>> th.runAnalysis(velocity, density, pSetting)
        >>> th.printResults()
        Thrust       = 1124.8730
        Power        = 73496.2684
        fuelFlow     = 0.0054
        """
        assert len(self.engineList)==len(self.propellerList)
        T      =np.zeros([self.nEng])
        P      =np.zeros([self.nEng])
        FF     =np.zeros([self.nEng])
        RPM    =np.zeros([self.nEng])
        N      =np.zeros([self.nEng])
        B      =np.zeros([self.nEng])
        for i,engine in enumerate(self.engineList):
            fixedPitch,fixedBeta=self._checkBetaRange(self.propellerList[i])
            propeller=self.propellerList[i]
            Teng=thrustPerEngine(engine,propeller)
            data=self._getThrust(Teng,V,rho,powerSetting,fixedPitch,fixedBeta)
            T[i]  =data[0]
            P[i]  =data[1]
            FF[i] =data[2]
            RPM[i]=data[3]
            N[i]  =data[4]
            B[i]  =data[5]
        self.thrust            =np.sum(T)
        self.power             =np.sum(P)
        self.fuelFlow          =np.sum(FF)
        self.RPM               =RPM
        self.propEfficiency    =N
        self.beta              =B
        self.powerSetting      =np.ones(len(self.engineList))*powerSetting
        self.SAR_km_per_kg     =V/self.fuelFlow/1000.
    def matchRequiredThrust(self,V,rho,requiredThrust):
        """
        Finds the power setting, RPM, fuel flow of the engine(s) to match
        an input required thrust at given velocity and density.
        
        Parameters
        ----------
        
        V : float, m/s
            velocity
        rho : float, kg/m**3
            air density
        requiredThrust : float, N
            thrust required

        Notes
        -----
        
        Assumes symmetric thrust for multi-engine configurations
        
        Examples
        --------
        
        Single engine variable pitch - Thrust matching
    
        >>> velocity = 40.0
        >>> requiredThrust = 600.0
        >>> th.matchRequiredThrust(velocity, density, requiredThrust)
        >>> th.printResults()
        Thrust       = 600.0000
        Power        = 26376.4080
        fuelFlow     = 0.0015
        """
        T      =np.zeros([self.nEng])
        P      =np.zeros([self.nEng])
        FF     =np.zeros([self.nEng])
        RPM    =np.zeros([self.nEng])
        N      =np.zeros([self.nEng])
        B      =np.zeros([self.nEng])
        PS     =np.zeros([self.nEng])
        assert len(self.engineList)==len(self.propellerList)        
        reqThrustPerEng=requiredThrust/len(self.engineList) # symmetric thrust

        for i,engine in enumerate(self.engineList):
            fixedPitch,fixedBeta=self._checkBetaRange(self.propellerList[i])
            propeller=self.propellerList[i]
            Teng=thrustPerEngine(engine,propeller)
            data=self._matchThrust(Teng,V,rho,reqThrustPerEng,fixedPitch,fixedBeta)
            PS    =data[0]
            T[i]  =data[1]
            P[i]  =data[2]
            FF[i] =data[3]
            RPM[i]=data[4]
            N[i]  =data[5]
            B[i]  =data[6]
        self.powerSetting      =PS     
        self.thrust            =np.sum(T)
        self.power             =np.sum(P)
        self.fuelFlow          =np.sum(FF)
        self.RPM               =RPM
        self.propEfficiency    =N
        self.beta              =B
        self.SAR_km_per_kg     =V/self.fuelFlow/1000.
    def _zeroThrustModule(self,V,rho,powerSetting):
        self.thrust            =0.
        self.power             =0.
        self.fuelFlow          =0.
        self.RPM               =0.
        self.propEfficiency    =0.
        self.beta              =0.
        self.powerSetting      =0.
        self.SAR_km_per_kg     =0. 
    def printResults(self):
        """
        Prints out thrust, power and fuel flow.
        """
        print '{0:12} = {1:.4f}'.format('Thrust',self.thrust)
        print '{0:12} = {1:.4f}'.format('Power',self.power)
        print '{0:12} = {1:.4f}'.format('fuelFlow',self.fuelFlow)
    def _checkBetaRange(self,propeller):
        """
        Checks propeller data from db.
        """
        if propeller.betaRange[0]==propeller.betaRange[1]:
            fixedPitch=True
            fixedBeta=propeller.betaRange[0]
        elif abs(propeller.betaRange[0]-propeller.betaRange[1])<=10:
            print 'Beta range is too low - setting fixed pitch'
            fixedBeta=.5*(propeller.betaRange[0]+propeller.betaRange[1])
            fixedPitch=True
        else:
            fixedPitch=False
            fixedBeta=None
        return fixedPitch,fixedBeta
    def _matchThrust(self,Teng,V,rho,requiredThrust,fixedPitch,fixedBeta):
        if fixedPitch==True:     
            data=Teng.thrustMatch_fixedPitch(V,rho,fixedBeta,requiredThrust)
        elif fixedPitch==False:
            data=Teng.thrustMatch_variablePitch(V,rho,requiredThrust)
        else:
            print 'attribute fixedPitch must be True or False'
            data=None
        return data
    def _getThrust(self,Teng,V,rho,powerSetting,fixedPitch,fixedBeta):
        if fixedPitch==True:     
            data=Teng.thrust_fixedPitch(V,rho,powerSetting,fixedBeta)
        elif fixedPitch==False:
            data=Teng.thrust_variablePitch(V,rho,powerSetting,mode=self.mode)
        else:
            print 'attribute fixedPitch must be True or False'
            data=None
        return data

#-------------THRUST MATCHING OBJECTIVE AND CONSTRAINT FUNCTIONS-------------
def _thrustMatchFixedObj(x,thrustObj,V,rho,beta,reqThrust,thrustFcn):
    RPM=x[0]
    powerSetting=x[1]
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho)   
    deltaT=(thrust-reqThrust)**2
    return deltaT
def _thrustMatchFixedCon(x,thrustObj,V,rho,beta,reqThrust,thrustFcn):
    RPM=x[0]
    powerSetting=x[1]
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho)   
    return deltaP
def _thrustMatchVariableObj(x,thrustObj,V,rho,reqThrust,thrustFcn):
    RPM=x[0]
    beta=x[1]
    powerSetting=x[2]
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho)       
    return -effy
def _thrustMatchVariableCon1(x,thrustObj,V,rho,reqThrust,thrustFcn):
    # Preq = Pavail
    RPM=x[0]
    beta=x[1]
    powerSetting=x[2]
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho)   
    return deltaP
def _thrustMatchVariableCon2(x,thrustObj,V,rho,reqThrust,thrustFcn):
    # Treq = Tavail
    RPM=x[0]
    beta=x[1]
    powerSetting=x[2]
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho) 
    deltaT=thrust-reqThrust
    return deltaT
#-------------THRUST AVAILABLE OBJECTIVE AND CONSTRAINT FUNCTIONS------------                
def _fixedPitchObj(x,thrustObj,V,rho,powerSetting,beta,thrustFcn):
    RPM=thrustObj.engine._checkRPM(x)
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho)    
    return deltaP**2

def _varPitchObj_maxEffy(x,thrustObj,V,rho,powerSetting,thrustFcn):
    RPM=thrustObj.engine._checkRPM(x[0])
    beta=x[1]
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho) 
    return -effy

def _varPitchObj_maxThrust(x,thrustObj,V,rho,powerSetting,thrustFcn):
    RPM=thrustObj.engine._checkRPM(x[0])
    beta=x[1]
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho)
    return -thrust
    
def _varPitchCon(x,thrustObj,V,rho,powerSetting,thrustFcn):
    RPM=thrustObj.engine._checkRPM(x[0])
    beta=x[1]
    thrust,power,deltaP,propRPM,effy=thrustFcn(thrustObj,RPM,beta,powerSetting,V,rho) 
    return deltaP
#----------------------------------------------------------------------------
def _getThrustData_exact(thrustObj,RPM,beta,throttle,V,rho):
    propRPM  =RPM/thrustObj.engine.gearboxRatio
    data     =thrustObj.propeller.analyze_prop(beta,propRPM,V,rho)
    if data[0]==None or np.isnan(data[0])==True:
        thrust   =-1.
        power    =-1.
        deltaP   =-1.
        effy     =-1.
    else:
        Pav      =thrustObj.engine.get_availablePower(rho,RPM,throttle)
        Pr       =data[1]
        deltaP   =Pav-Pr
        thrust   =data[0]
        effy     =data[5]
        power    =Pr
    return thrust,power,deltaP,propRPM,effy

def _getThrustData_approx(thrustObj,RPM,beta,throttle,V,rho):
    propRPM  =RPM/thrustObj.engine.gearboxRatio
    n      =propRPM/60.
    d      =thrustObj.propeller.diameter
    Pa     =thrustObj.engine.get_availablePower(rho,RPM,throttle)
    CP     =thrustObj.propTable.getCP(beta,propRPM,V)
    CT     =thrustObj.propTable.getCT(beta,propRPM,V)
    effy   =thrustObj.propTable.getEffy(CT,CP,propRPM,V)
    thrust =CT*rho*n**2*d**4
    Pr     =CP*rho*n**3*d**5    
    deltaP =Pa-Pr
    power  =Pr
    return thrust,power,deltaP,propRPM,effy
class thrustPerEngine():
    def __init__(self,engine,propeller,propTable=None):
        self.engine      =engine
        self.propeller   =propeller
        self.propTable   =propTable
    def _selectThrustFunction(self,forceExactMode):
        if self.propTable==None or forceExactMode==True:
            thrustFun=_getThrustData_exact
        else:
            thrustFun=_getThrustData_approx
        return thrustFun
    def _checkResults(self,residual,data,tolerance):
        thrust        =data[0]
        power         =data[1]
        deltaP        =data[2]
        RPM           =data[3]
        effy          =data[4]
        if residual <-tolerance or thrust<0 or np.isnan(thrust) or np.isnan(power):
            #print 'Warning - no converged solution'
            thrust    =None
            power     =None
            deltaP    =None
            RPM       =None
            effy      =None
        return thrust,power,deltaP,RPM,effy
    def setPropTable(self,propTable):
        self.propTable =propTable
    def buildPropTable(self,Blist,propRPMlist,Vlist):
        self.propTable=propTable(self.engine,self.propeller,Blist,propRPMlist,Vlist)
    def thrust_fixedPitch(self,V,rho,powerSetting,beta,forceExactMode=True):
        RPMLo            =self.engine.get_minRPM()
        RPMHi            =self.engine.get_maxRPM()
        fun              =_fixedPitchObj
        bnds             =[[RPMLo,RPMHi]]
        thrustFun        =self._selectThrustFunction(forceExactMode)
        Args             =self,V,rho,powerSetting,beta,thrustFun
        tolerance        =10**-3       
        residual         =999.0
        result           =minimize_scalar(fun, bracket=bnds[0], bounds=bnds[0], args=Args, method='Bounded')
        residual         =result.fun
        RPM              =self.engine._checkRPM(result.x)
        data             =thrustFun(self,RPM,beta,powerSetting,V,rho)
        thrust,power,deltaP,propRPM,effy=self._checkResults(residual,data,tolerance)
        fuelFlow         =self.engine.get_fuelFlow(rho,RPM,powerSetting)
        return thrust,power,fuelFlow,RPM,effy,beta
    def thrustMatch_fixedPitch(self,V,rho,beta,thrustReq,forceExactMode=True):
        RPMLo            =self.engine.get_minRPM()
        RPMHi            =self.engine.get_maxRPM()
        fun              =_thrustMatchFixedObj
        bnds             =[[RPMLo,RPMHi],[10.,100.]]
        x0               =[(RPMLo+RPMHi)/2.,50.]
        thrustFun        =self._selectThrustFunction(forceExactMode)
        args             =self,V,rho,beta,thrustReq,thrustFun
        cons             =({'type':'eq','fun':_thrustMatchFixedCon,'args':args})
        opts             ={'maxiter':100,'disp':False}
        tolerance        =10**-3
        residual         =999.0
        result           =minimize(fun,x0,args,method='SLSQP',bounds=bnds,constraints=cons,options=opts)
        residual         =result.fun
        RPM              =self.engine._checkRPM(result.x[0])
        powerSetting     =result.x[1]
        data             =thrustFun(self,RPM,beta,powerSetting,V,rho)
        thrust,power,deltaP,propRPM,effy=self._checkResults(residual,data,tolerance)
        fuelFlow         =self.engine.get_fuelFlow(rho,RPM,powerSetting)
        return powerSetting,thrust,power,fuelFlow,RPM,effy,beta
    def thrust_variablePitch(self,V,rho,powerSetting,mode='max_thrust',forceExactMode='True'):
        if mode=='max_thrust':
            fun=_varPitchObj_maxThrust
        elif mode=='max_efficiency':
            fun=_varPitchObj_maxEffy       
        else:
            fun=_varPitchObj_maxThrust
            print 'invalid mode - reverting to maximum thrust mode'
            print 'please specify either max_thrust or max_efficiency'
        betaLo           =self.propeller.betaRange[0]
        betaHi           =self.propeller.betaRange[1]
        RPMLo            =self.engine.get_minRPM()
        RPMHi            =self.engine.get_maxRPM()       
        bnds             =[[RPMLo, RPMHi], [betaLo, betaHi]]
        x0               =[(RPMLo+RPMHi)/2.,(betaLo+betaHi)/2.]
        thrustFun        =self._selectThrustFunction(forceExactMode)
        args             =self,V,rho,powerSetting,thrustFun
        cons             =({'type':'eq','fun':_varPitchCon,'args':args})
        opts             ={'maxiter':100,'disp':False}
        tolerance        =10**-3                
        residual         =0.0  
        res              =minimize(fun,x0,args,method='SLSQP',bounds=bnds,constraints=cons,options=opts)
        RPM              =res.x[0]
        beta             =res.x[1]
        data             =thrustFun(self,RPM,beta,powerSetting,V,rho)
        thrust,power,deltaP,propRPM,effy=self._checkResults(residual,data,tolerance)
        fuelFlow         =self.engine.get_fuelFlow(rho,RPM,powerSetting)
        return thrust,power,fuelFlow,RPM,effy,beta
    def thrustMatch_variablePitch(self,V,rho,thrustReq,forceExactMode=True):
        betaLo           =self.propeller.betaRange[0]
        betaHi           =self.propeller.betaRange[1]
        RPMLo            =self.engine.get_minRPM()
        RPMHi            =self.engine.get_maxRPM()
        fun              =_thrustMatchVariableObj
        bnds             =[[RPMLo,RPMHi],[betaLo,betaHi],[10.,100.]]
        x0               =[(RPMLo+RPMHi)/2.,(betaLo+betaHi)/2.,50.]
        thrustFun        =self._selectThrustFunction(forceExactMode)
        args             =self,V,rho,thrustReq,thrustFun
        cons             =({'type':'eq','fun':_thrustMatchVariableCon1,'args':args},{'type':'eq','fun':_thrustMatchVariableCon2,'args':args})
        opts             ={'maxiter':100,'disp':False}
        tolerance        =10**-3                
        residual         =999.0               
        result           =minimize(fun,x0,args,method='SLSQP',bounds=bnds,constraints=cons,options=opts)
        residual         =0.0
        RPM              =self.engine._checkRPM(result.x[0])
        beta             =result.x[1]
        powerSetting     =result.x[2]
        data             =thrustFun(self,RPM,beta,powerSetting,V,rho)
        thrust,power,deltaP,propRPM,effy=self._checkResults(residual,data,tolerance)
        fuelFlow         =self.engine.get_fuelFlow(rho,RPM,powerSetting)
        return powerSetting,thrust,power,fuelFlow,RPM,effy,beta
class propTable():
    def __init__(self,engine,propeller,Blist,propRPMlist,Vlist):
        atm=ISAtmosphere(0.0)
        self.engine      =engine
        self.propeller   =propeller
        self.rhoSSL      =atm.density
        self.Blist       =Blist
        self.propRPMlist =propRPMlist
        self.Vlist       =Vlist
        self.buildTalbes()
    def buildTalbes(self):
        mesh_B=np.zeros([len(self.Blist)*len(self.propRPMlist)*len(self.Vlist)])
        mesh_R=np.zeros([len(self.Blist)*len(self.propRPMlist)*len(self.Vlist)])
        mesh_V=np.zeros([len(self.Blist)*len(self.propRPMlist)*len(self.Vlist)])
        mesh_P=np.zeros([len(self.Blist)*len(self.propRPMlist)*len(self.Vlist)])
        mesh_T=np.zeros([len(self.Blist)*len(self.propRPMlist)*len(self.Vlist)])
        q=0
        for i,V in enumerate(self.Blist):
            for j,RPM in enumerate(self.propRPMlist):
                for k,B in enumerate(self.Vlist):
                     CT,CP,runFail=self._getResults(B,RPM,V)  
                     if runFail==False:
                         mesh_B[q]=B
                         mesh_R[q]=RPM
                         mesh_V[q]=V
                         mesh_P[q]=CP
                         mesh_T[q]=CT
                         q+=1
        self.mesh_B=mesh_B[0:q-1]
        self.mesh_R=mesh_R[0:q-1]
        self.mesh_V=mesh_V[0:q-1]
        self.mesh_P=mesh_P[0:q-1]
        self.mesh_T=mesh_T[0:q-1]
        self.points=np.zeros([len(self.mesh_B),3]) 
        self.points[:,0]=self.mesh_B
        self.points[:,1]=self.mesh_R
        self.points[:,2]=self.mesh_V
        self.interpCP=interp.LinearNDInterpolator(self.points,self.mesh_P)
        self.interpCT=interp.LinearNDInterpolator(self.points,self.mesh_T)
    def _getResults(self,beta,propRPM,V):
        data=self.propeller.analyze_prop(beta,propRPM,V,1.225)
        CT=data[3]
        CP=data[4]
        runFail=False
        if data[0]==None:runFail=True
        return CT,CP,runFail
    def getCP(self,beta,propRPM,V):
        return self.interpCP(beta,propRPM,V)
    def getCT(self,beta,propRPM,V):
        return self.interpCT(beta,propRPM,V)  
    def getEffy(self,CT,CP,propRPM,V):
        n=propRPM/60.0
        J=V/n/self.propeller.diameter
        N=J*CT/CP
        return N
    def getPropPower(self,beta,propRPM,V,rho):
        CP=self.getCP(beta,propRPM,V)
        n=propRPM/60.
        d=self.propeller.diameter
        return CP*rho*n**3*d**5
    def getPropThrust(self,beta,propRPM,V,rho):
        CT=self.getCT(beta,propRPM,V)
        n=propRPM/60.
        d=self.propeller.diameter
        return CT*rho*n**2*d**4


class totalThrust():
    def __init_(self,aircraft):
        pass
class asymmetricThrust():
    def __init__(self,aircraft,throttleSettingList):
        pass

#                 ---TEST SECTION---

def runTest_varPitch():
    import aircraft
    import matplotlib.pyplot as plt
    import matplotlib
    ac=aircraft.load("V0510")
    T=thrustPerEngine(ac.engine.engineList[0],ac.propeller.propellerList[0])
    V           =76.
    rho         =1.2
    pctThrottle =75.   
    
#    N=20
#    Blist       =frange(15.,45.,N)
#    propRPMlist =frange(500.,3000.,N)
#    Vlist       =frange(1,70,N)       
    #T.buildPropTable(Blist,propRPMlist,Vlist)   

    Vlist =np.array(range(30,70,2))
    thrust1=np.zeros(len(Vlist))
    effy1  =np.zeros(len(Vlist))
    beta1  =np.zeros(len(Vlist))
    thrust2=np.zeros(len(Vlist))
    effy2  =np.zeros(len(Vlist))
    beta2  =np.zeros(len(Vlist))
    
    t1 = time.time()
    for i in range(0,100,1):
        T.thrust_variablePitch(V,rho,pctThrottle,mode='max_thrust')
    t2 = time.time()
    print 'Time = '+str((t2-t1)/100)+ ' seconds per run'
    
    for i,V in enumerate(Vlist):
        thrust1[i],powerAchieved1,fuelFlow1,RPM,effy1[i],beta1[i]=T.thrust_variablePitch(V,rho,pctThrottle,mode='max_thrust')
        thrust2[i],powerAchieved2,fuelFlow2,RPM,effy2[i],beta2[i]=T.thrust_variablePitch(V,rho,pctThrottle,mode='max_efficiency')
        print str(thrust1[i])+' '+str(thrust2[i])
    
    font = {'family' : 'sans-serif',
            'style'  : 'normal',
            'size'   : 16}
    matplotlib.rc('font', **font)
    
    plt.figure(1)
    plt.hold(True)
    plt.grid(True)
    l1,=plt.plot(Vlist,effy1,'b-')    
    l2,=plt.plot(Vlist,effy2,'r-')
    plt.legend([l1,l2],['Max Thrust Mode','Max Efficiency Mode'],loc=2)
    plt.title('Velocity vs. Efficiency')
    plt.xlabel('Velocity (m/s)')
    plt.ylabel('Efficiency')
    
    plt.figure(2)
    plt.hold(True)
    plt.grid(True)
    l1,=plt.plot(Vlist,thrust1,'b-')
    l2,=plt.plot(Vlist,thrust2,'r-')
    plt.legend([l1,l2],['Max Thrust Mode','Max Efficiency Mode'],loc=2)
    plt.title('Velocity vs. Thrust')
    plt.xlabel('Velocity (m/s)')
    plt.ylabel('Thrust (N)')
    
    plt.figure(3)
    plt.hold(True)
    plt.grid(True)
    l1,=plt.plot(Vlist,beta1,'b-')
    l2,=plt.plot(Vlist,beta2,'r-')
    plt.legend([l1,l2],['Max Thrust Mode','Max Efficiency Mode'],loc=2)
    plt.title('Velocity vs. Beta')
    plt.xlabel('Velocity (m/s)')
    plt.ylabel('Beta (deg)')

    plt.show()
def runTest_fixPitch():
    import aircraft
    import matplotlib.pyplot as plt 
    ac          =aircraft.load("sampleInput1")
    thrust      =thrustPerEngine(ac.engine.engineList[0],ac.propeller.propellerList[0])
    
    N=10
    Blist       =frange(15.,45.,N) # pitch angle list
    propRPMlist =frange(500.,3000.,N)
    Vlist       =frange(1,70,N)       
    thrust.buildPropTable(Blist,propRPMlist,Vlist)
    
    V           =30.
    rho         =1.2
    beta        =20.0

    #pwrList =np.array(frange(50.,100.,20))
    pwrList = np.array([50.,75.,100.])
    Vlist   =np.array(frange(1., 70.,20))

    t1 = time.time()
    for i in range(0,100,1):thrust.thrust_fixedPitch(V,rho,100,beta,True)
    t2=time.time()    
    for i in range(0,100,1):thrust.thrust_fixedPitch(V,rho,100,beta,False)
    t3 = time.time()
    print 'exact  time = '+str(t2-t1)
    print 'approx time = '+str(t3-t2)
    
    T  =np.zeros([len(pwrList),len(Vlist)])
    P  =np.zeros([len(pwrList),len(Vlist)])
    N  =np.zeros([len(pwrList),len(Vlist)])
    RPM=np.zeros([len(pwrList),len(Vlist)])
    VV =np.zeros([len(pwrList),len(Vlist)])
    
    beta  =25.    
    plt.figure(1)
    plt.hold(True)
    for i,pwr in enumerate(pwrList):
        for j,V in enumerate(Vlist):
            T[i,j],P[i,j],FF,RPM[i,j],N[i,j],beta=thrust.thrust_fixedPitch(V,rho,pwr,beta,1)
            #aprx=thrust.thrust_fixedPitch(V,rho,pwr,beta,0)            
            VV[i,j]=V
        plt.plot(VV[i,:],T[i,:],'o-')
    plt.xlabel('Velocity, m/sec')
    plt.ylabel('Thrust, N')
    plt.grid(True)
    plt.legend(['50% throttle','75% throttle','full throttle'])
    plt.show()
def runTest_propTable():
    import aircraft
    import matplotlib.pyplot as plt 
    N=20
    ac          =aircraft.load("sampleInput1")
    prop        =ac.propeller.propellerList[0]
    Blist       =frange(15.,45.,N)
    propRPMlist =frange(500.,3000.,N)
    Vlist       =frange(1,70,N)
        
    thrust      =thrustPerEngine(ac.engine.engineList[0],ac.propeller.propellerList[0])
    thrust.buildPropTable(Blist,propRPMlist,Vlist)            
    exactMode =thrust.thrust_fixedPitch(40.,1.22,100,25.,True)  
    approxMode=thrust.thrust_fixedPitch(40.,1.22,100,25.,False)  
    print exactMode
    print approxMode
    #propTbl     =propTable(engine,prop,Blist,propRPMlist,Vlist)
    propTbl      =thrust.propTable
    
    beta    =25.
    rho     =1.225
    propRPM =1500.
    V       =35.
    print 'CP     = '+str(propTbl.getCP(beta,propRPM,V))
    print 'Power  = '+str(propTbl.getPropPower(beta,propRPM,V,rho)/1000.)+' kw'
    print 'Thrust = '+str(propTbl.getPropThrust(beta,propRPM,V,rho))+' N'
    
    T,P,Jreq,CT,CP,effy=prop.analyze_prop(beta,propRPM,V,1.225)
    print 'CP     = '+str(CP)
    print 'Power  = '+str(P/1000.)+' kw'
    print 'Thrust = '+str(T)+' N'
    
    T=np.zeros([len(Vlist)])

    plt.figure(1)
    plt.hold(True)
    for i,RPM in enumerate(propRPMlist):
        for j,V in enumerate(Vlist):
            #T[j]=CPtab.getPropThrust(beta,RPM,V,rho)
            T[j]=propTbl.getCT(beta,RPM,V)
        plt.plot(Vlist,T,'o-')
    plt.show()
def runTest_thrustAnalysis():
    import aircraft
    print '\n Testing Single Engine Variable Pitch'
    ac=aircraft.load("sampleInput1")
    th=thrustAnalysis(ac)
    th.runAnalysis(60.,1.225,100.)
    th.printResults()
    print '\n Testing Thrust Matching - Single Engine Variable Pitch'
    th.matchRequiredThrust(40.,1.225,600)
    th.printResults()
    print '\n Testing - Twin Engine Variable Pitch'
    ac=aircraft.load("sampleInput2")
    th=thrustAnalysis(ac)
    th.runAnalysis(60.,1.225,100.)
    th.printResults()
    print '\n Testing Thrust Matching - Twin Engine Variable Pitch'
    th.matchRequiredThrust(40.,1.225,1200)
    th.printResults()
    print '\n Testing Twin Engine Fixed Pitch'
    ac=aircraft.load("fixedPitchSample")
    th=thrustAnalysis(ac)
    th.runAnalysis(60.,1.225,100.)
    th.printResults()
    print '\n Testing Thrust Matching - Twin Engine Fixed Pitch'
    th.matchRequiredThrust(40.,1.225,1200)
    th.printResults()
if __name__=="__main__":
    runTest_thrustAnalysis()
    runTest_varPitch()
    runTest_fixPitch()
    #runTest_varPitch()
    #runTest_fixPitch()
    #runTest_propTable()
 